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Abstract. We study the two-qubit Rabi model in the most general case where the 
qubits are different from each other. The spectrum of the system in the ultrastrong- 
coupling regime is shown to converge to two forced oscillator chains by perturbation 
theory. An even and odd decomposition of the Hilbcrt space allows us to calculate 
the spectra in any given parameter regime; the cases studied confirm our perturbation 
theory prediction in the ultrastrong-coupling regime and point to crossings in the 
spectra within each parity subspacc in the moderate-coupling regime. The normal 
modes of the system are calculated by two different methods, the first a linear algebra 
approach via the parity bases that delivers a four-term recurrence relation for the 
amplitudes of the proper states and, the second, via Bargmann representation for the 
field that delivers five-term recurrence relations. Finally, we show some examples of 
the time evolution of the mean photon number, population inversion, von Neuman 
entropy and Wootters concurrence under the two-qubit quantum Rabi Hamiltonian by 
taking advantage of the parity decomposition. 
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1. Introduction 

The Dicke model [I] is the simplest model describing the interaction of an ensemble of 
identical two-level systems with a boson field under the rotating wave approximation. 
In units of H, it is written as: 

Hd = Ufa^a + ojqJz + g (^aJ^ + a}j_^ , (1) 

where the field mode of frequency is described by creation (annihilation) operators 
(a) and the identical two-level systems or qubits with transition frequency uq are 
described by collective Dicke or angular momentum operators answering the su{2) 
algebra, [Jz,j±] = ±</± and [J+, J_] = 27^, the ensemble-field coupling is given by 
the parameter g. The Dicke model, also known as Tavis-Cummings model [2], is exactly 
solvable and has successfully used to describe super-radiance, super-fiuorescence and 
amplified spontaneous emission phenomena in cavity-quantum electrodynamics (cavity- 
QED), ion-trap cavity-QED, and circuit-QED systems; see [3] for a review. 
Without the rotating wave approximation (RWA) and in units of h, 

Hr = UfO^a + UqJz + g [a + a^) + J_ j , (2) 

the Dicke model has recently been known as the quantum Rabi model. For many years 
the quantum Rabi model was a mathematical curiosity, as dipole-field systems cannot 
reach the strong-coupling regime in a simple way; e.g., the proposed and experimental 
realizations to find the quantum phase transition in the ground state of the Dicke model 
[H El [HI [7] . Nevertheless, the proper system, dynamics and the notion of integrability of 
a single two-level system interacting with a single mode boson field without the RWA 
was studied utihzing Bargman representation [8] , continued fractions [9l [TOl [11] , semi- 
classical methods [12], perturbation theory [131 E], and coupled cluster [TS] methods. It 
has also been discussed in the literature the influence of the so-called counter-rotating 
terms eliminated by the RWA on the phase distribution of a micromaser cavity field 
[16j . on the quantum Zeno and anti-Zeno effects [T^ and on the linear polarizability 
[I8] . Recently, the interest on the single qubit Rabi model has been relighted due to 
the ability of solid-state-cavity-QED P^] and circuit-QED [201 121] systems to attain 
the strong-coupling regime where the RWA fails. Thus, the single qubit Rabi model has 
been revisited using perturbation theory [221123], continued fractions [211 125], Bargmann 
representation [261 EH EHl EHl [30l [31] , and coherent states [32] methods. 

Here, we are interested in the quantum Rabi model for two nonidentical qubits, 
described by the Pauli operators a^p with j = x,y,z and = 1,2 and the transition 
frequencies ui and U2, 

H = cofd^a + ^ {co.aW + u:,aP) + {d + at) {g,a^^) + g,a^^)) , (3) 

which, up to our knowledge, has only been studied in two simplified forms: one for 
identical qubits with frequencies much smaller than the oscillator frequency [33] and 
the other for non-identical qubits where the transition frequency of one of them is set 
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to zero [S]. The model Q conserves parity defined as 

fI = a«af(-l)'^^^ (4) 

1. e., [n, H] = 0, and may be of interest to tlie circuit-QED community wliere it lias been 
used in the weak-coupling regime under the RWA for coherent quantum storage and 
transfer between two phase qubits [35] and for resonant two-qubit phase gates the 
validity of the latter breaks in the ultrastrong-coupling regime where the full quantum 
Rabi model has to be used |37|. The effect of counter rotating terms on entanglement 
and discord in a lossy two-qubit Rabi system has also been studied numerically [55] . 

In the following, we discuss the eigenvalues of Hamiltonian ^ in both the weak- and 
strong-coupling regimes. The eigenvalues are exact in the former and approximated up 
to second order correction by using perturbation theory in the latter. We take advantage 
of parity conservation to numerical calculate the spectra in each of the parity subspaces. 
Then, we present the proper functions via linear algebra methods and functions in 
Bargmann representation [39]. In the second to last section, we show how the time 
evolution in the weak-coupling regime can be calculated in closed form and how to 
calculate the evolution of quantities of interest in the parity decomposition. Finally, we 
close with a brief conclusion. 

2. Eigenvalues of the model 

In this section we calculate the exact eigenvalues for the system in the weak-coupling 
regime by partitioning the Hilbert space into subspaces described by the total excitation 
number. In the strong-coupling regime, we use perturbation theory to approximate the 
eigenvalues to second order correction. 

2.1. Weak-coupling regime 

It is quite trivial to calculate the eigenvalues of the general two-qubit quantum Rabi 
model ([3]) in the weak-coupling regime, gi^g2 ^ ujf, where the RWA approximation is 
valid, 

i=l,2 j=l,2 

In this regime, we can find the eigenvalues by partitioning the Hilbert space into 
subspaces with a well defined total excitation number, N = n + (^ai^^ + a^^^^ /2; for 
example, for zero excitation the eigenstate is |0, g, g) with energy Aq = — (wi + u;2)/2, for 
one excitation the subspace is span by {|1, g, g), |0, g, e), |0, e, g)}, and for the rest, n > 2, 
the corresponding subspace is span by {\n, g, g),\n — 1, g, e),\n — 1, e, g),\n — 2, e, e)} 

Figure ([T]) shows the exact spectra in the RWA versus the numerical spectra for 
the whole Hamiltonian considering up to one thousand excitations in the system. We 
present in Fig. [T]^a) the on-resonance case for equal qubits, coi = 102 = ojf, and equal 
couplings gi = g2 = g G [0,0.1510;/. Figure [T]^b) shows the case presented in [33], where 
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one of the qubits has zero self-energy and the quantum field is detuned to the blue of the 
qubit with nonzero transition frequency, Ui = O.Tuf and U2 = 0, and equal couplings, 
91 = 92 = 9 ^ [0, O.lSjcj/. The case for symmetric detuning, ui = 0.9uf and CO2 = l.lcof, 
with the first qubit coupled to the field at a fixed coupling value gi = O.Olujf and variable 
coupling parameter for the other qubit, 92 = 9 & [0,0.15]^/, is shown in FigQc). We 
show a similar case, ui = l-lujf and 002 = 0.9a;/, (72 = 9 ^ [0, 0.15]a;/, with a slightly 
larger fixed couphng, gi = O.lcof in Fig. ^d). One can see that it is safe to say that 
the exact eigenvalues in the RWA describe well the spectra when the relation between 
the coupling strengths and the field frequency is below 10%, 9j/ujf < 0.1, [37] but some 
particular cases may extend further this range in one coupling by restricting the other; 
e.g.. Figure [l](c) and[l|^d). Notice that the RWA fails to describe the anti-crossings in 
the quantum Rabi spectra; e.g.. Figure [l|^a) and[l](d). 



2.2. Strong- coupling regime 

In order to approximate the eigenvalues in this regime, let us assume that the coupling 
parameters are larger than the transition frequencies, gi,g2 ^ (^1,(^2- In such a case, 
we can write the two-qubit quantum Rabi model (|3]) as a leading Hamiltonian with a 
perturbartion, H = Hq + P, where 



Ho = Uffi + (a + a^) {giaj^^^ + g2a^^^) 
At this point we can use the transformation 



(6) 
(7) 



— ^ - Yy ^ - (8) 

with 9 = Tt/A to implement the changes ai"''' — )■ ai-'^ and cri"''' — ?■ —ctx^ leading to the 
form, 



^0 = ^ffi + (a + a^) {gia^^^ + (72^^)) , 



(9) 

(10) 



where the tilde has been used to represent the rotated operator, O = K^y{Tc / 4)0 Ry{TT / A) . 
Notice that the unperturbed Hamiltonian Ho can be diagonalized, 

\ 



H^ = T, 



D 



( 9l 



n — 



\ 



^2 
^2 
0^-2 



:iii 



D 



9iJ 

by using a driven oscillator basis provided by the transformation 

/ D{g+lujf) 

D{g_/uf) 

D{-g^/ujf) 

V D{-9+/uJf) J 



(12) 




Figure 1. (Color online) Exact spectra in the rotating- wave approximation (dashed 
blue) and numerical spectra for the quantum Rabi Hamiltonian considering a 
subspace of up to one thousand excitations (solid red) for different parameter sets: 
91.92} = (a) {l,l,5,g}L./, (b) {0.7, 0, 5, gjc./, (c) {0.9, 1.1, 0.01, g^, and (d) 
{l.l,0.9,0.1,g}w/. 

where the displacement operator is defined as D{a) = e°'"-''~°''°- fulfilhng /)"''(«) = D{—a) 
and D{a)D{/3) = D{a + /3). The auxihary couphng parameters are defined as 
g± = gi ± g2. Thus, we have the zeroth order approximated eigenvalues in the strong- 
coupling regime, 

^2 

ei,m = e4,m ~ w/m (13) 

Uf 

e2,m = e3,m ~ - — , (14) 
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that happen to be two-fold degenerate. The first order correction due to the perturbation 
P is null due to the lack of diagonal components in the driven oscillator basis, 

/ U2D{2g2/ujf) uJiD{2gi/uf) 

^ ^ i02D{-2g2/ujf) uJiD{2g,/iOf) 

^ ujrb{-2g^/ujf) U2b{2g2/ujf) 

\ uJib{-2gi/ujf) uj2b{-2g2/ujf) 



but the second order corrections are available leading to the approximate eigenvalues, 

ei,m = e4,m ~ ojfm - ^ - e2,+, (16) 

Uf 

e2,m = e3,m ~ c^/"^ - — - £2,-, (17) 

Uf 

where the second order correction is given by 

,^^=y u^!\{m\b{2gi/uf)\n)\' ^ y uj\{m\bi2g2/ujf)\n)\' 
^'^ ^^/("^-^) ±%i^2M ^^UJf{m-n)±4gig2/ujf' 

Notice that for extremely large values of the coupling parameters gi and g2, i.e., 
ultrastrong-coupling, the second order corrections tend to zero and each of the 
spectral branches of the two-qubit quantum Rabi model, (16) and (17), tends to 
the spectra of a forced oscillator. For large values of the coupling parameters 
91, 92 ^ ^1, 1^2, we can calculate the second order correction by using the identity 
{m\b{2x)\n) = -y^(2x)("-")e-2|^l'£i™"''^(|2xn where C'a\z) is the 6th generalized 
Laguerre polynomial. 

Figure [2] shows the numerical spectra for the quantum Rabi Hamiltonian expanded 
in the parity subspaces span by the even, (21), and odd, (22), parity basis defined in 
the following section. We consider up to a thousand excitations in the system; i.e., 
truncated H± matrices given by (23) of size one thousand with coupling steps of one 
hundredth of the field frequency, Ag = 0.01 cOf. In Figj2|^a) and (b), we show the 
symmetric off-resonance case, Ui = 1.3 uj and U2 = 0.7 uj, with identical couplings, 
9i = 92 = 9 ^ [0)2]c<;/. The spectra corresponding to even parity eigenvectors is 
presented as solid red lines and that corresponding to odd parity as dashed blue lines. 
The inset in Fig. |2](a) shows an energy crossing in the even spectra confirmed by 
comparison of the proper states before and after the crossing. The spectra in the 
moderate coupling region is presented in Fig. [2|b). 

Figure [2]^c) and (d) shows the numerical spectra for the quantum Rabi Hamiltonian 
for the case presented in [31]; i-e., blue detuned field with the transition frequency 
of one qubit being negligible, ui = 0.7a;/ and ijJ2 = 0, and identical couplings, 
9i = 92 = 9 ^ ['^,2]ujf. Again, sohd red lines correspond to the spectra of the even 
parity subspace and dashed blue lines to that in the odd parity subspace. We present a 
magnification of this case in the moderate coupling region in|2](d). Note that the even 
and odd spectra are completely degenerate. Apparent energy crossings appeared in all 
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the cases studied but none was confirmed neither by decreasing the step size nor by 
comparing the eigenvectors before and after the apparent crossing; i.e., we looked at the 
first fifty eigenvalues by using a random transition frequency ui G {0,2]u)f. The inset 
in Fig. ^c) shows a typical energy anti-crossing in this particular case. 

We also studied a variety of cases not shown here; e.g. random qubit transition 
energies uji,ijJ2 G [0,2]ci;/ with identical couplings gi = g2 = g or one fixed random 
coupling gi G [0,2]a;/. In all the cases studied, energy crossings in the spectra of each 
parity subspaces were the norm and the spectra of the parity subspaces always cross 
each other. Also, in the cases of identical couplings, gi = g2 = g, the lower eigenvalue 
branches already tend to the approximated spectra of two driven oscillator chains for 
not-so-large values of the coupling parameters; in the cases shown in Fig. and ^c) 
the first dozen eigenvalue branches already show this behavior for gi = g2 > 2a;/. 



3. Eigenstates of the model 

Here we take on the calculation of the proper states of the system corresponding to the 
eigenvalues of the two-qubits quantum Rabi model. For the sake of completeness we use 
both a linear algebra and the Bargmann representation method that lead to four- and 
five-term recurrence relations between the coefficients of the eigenstates, respectively. 
The eigenstates provided by both methods are valid for any given coupling parameter 
value. 



3.1. Linear algebra approach 

In order to calculate the proper functions of the model (|3]), let us use the fact that it 
conserves parity defined in Q. Then, it is trivial to extend the idea of a parity basis 
used in the single-qubit Rabi model [TUl [TTl |22], 

\+,j) = {V^a^y\0)f\g)i, (19) 
|-,j) = (V^^^.)^|0)/|e)i, (20) 

where the Susskind-Glogower operator is defined as = YlT=o + PO]) such that 
we obtain a even and odd parity bases, 

= {|+,0,(?),|-,0,e),|+,l,^7),|-,l,e),|+,2,(?),|-,2,e),...},(21) 
= {|-,0,(7),|+,0,e),|-,l,(?),|+,l,e),|-,2,(7),|+,2,e),...},(22) 



in that order, with n|j)-|- = ±|j)±. The model Hamiltonian ([3j) in these bases, 
{H±)j^k =± {j\H\k)±, has block tridiagonal form: 

/ Oi . . . \ 
Oi Df O2 ... 
O2 O3 ... 



(23) 




Figure 2. (Color online) Numerical spectra for the quantum Rabi Hamiltonian 
considering subspaces of up to one thousand excitations for even (solid red) and odd 
(dashed blue) parity subspaces for different parameter sets: {wi, W2, 3i, 32} — (a) and 
(b) {1.3,0.7, g,g}ujf, (c) and (d) {0.7,0,g,g}ojf. The inset in (a) shows a crossing 
in the even parity spectra of the system at hand. The inset in (c) shows a typical 
anti-crossing in the corresponding configuration. The driven oscillator behavior of the 
spectra, predicted by perturbation theory for extremely large values of the coupling 
parameters with respect to the transition frequencies, can already be seen for couplings 
twice as large as the transition frequencies for the lower eigenvalues. 
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where the blocks are given by 



dt 




Df = I I , o, = I r I • (24) 

with d'^ = jujf =F I [(— ± UJ2] and = juf ± | [(— l)-'ci;i ± UJ2]. Notice that the 
matrices Oj are invertible as long as l^fip 7^ |5'2p- Thus, the eigenstates are written as 

00 

oc ^^W|+,j,^) + t;W|_,J,e), (25) 

j=0 

00 

le-) oc +4-)|+,j,e), (26) 

where we have used the symbol C,± to emphasize that the proper functions delivered by 
this method are valid for any given coupling parameter values. The coefficients can be 
calculated up to a normalization factor by the four-term recurrence relations given by 

4^)= - O^' {D^ - IQ vit\ (27) 
= - Oj' {Df_, - U^) v^^\ - 0.'0,^,^^i J =2,3,... (28) 
by choosing a suitable vl^\ with Vj^^ = {Vj^\vj^'') and using 

(29) 




where the symbol 1 is the unitary matrix of dimension two. 

Figure [3] shows the normalized ground state for the case of two nonidentical qubits 
off-resonance with the field frequency, Ui = 1.1 uj and U2 = 0.3 Uf, for two strong- 
coupling parameters, gi = 3 Uf and g2 = 4 Uf. For the sake of simplicity we have 
defined the normalized coefficients as 'U2j+fc °^ "^fk- ^^^^ '^^^^ ground state energy 
is degenerate for even and odd parity chains, ^+ = = —49.0052 huf. Numerical 
evidence points to the fact that for large coupling parameters, gi,g2 ^ oJj, the ground 
state energy is degenerate for even and odd parity chains independent of the qubit 
transition frequency values, ui and U2. We want to emphasize that keeping track of the 
relative error is of great importance when calculating the eigenstates via the recurrence 



relations given in (27) and (28). A safer approach is to utilize standard linear algebra 



diagonalization packages for large sparse matrices. 
3.2. Bargmann representation method 

We can follow the procedure presented in p4J with our Hamiltonian model (jsj). That 
is, we use the transformation ^ to arrive to the rotated Hamiltonian, 

Hr = RlHRy, (31) 
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Figure 3. (Color online) The (a) even and (b) odd parity ground state coefficients 
for the quantum Rabi Hamiltonian considering a subspace of up to one thousand 
excitations for the parameter set {a;i,a;2, 51,52} = {1.1, 0.3, 3, 4}a;/. The ground states 
are degenerate with energy ^± sa —49.0052 hojf. 



H 



R 



= tufd^a - ^ + L02a^^^) + {a + a^) {gi&i''^ + g2af^) , (32) 

In the basis defined by ai^^ and using the Bargmann representation method, i.e., making 
the substitution d ^ dz and d'^ ^ z where the shorthand notation dx has been used to 
denote partial derivations with respect to x, this rotated Hamiltonian looks like 

Ufzdz - f cr^) + ((?2ai'^ + gi){z + d,) -f 

-f ojfzdz - f ai') + ig^ai'^ - g,){z + d,) 

(33) 

Then, we can use the so-called Fulton-Gouterman transformation [lH |26l [31] , 

TpG ~ — 1= I (2) (2) ) 5 (34) 

V2 \ TpGCTx —TfgO'x J 

where the map Tpc is defined as Tpc : (p{z) — )■ 0(— 2;). This procedure yields 

flGHRTpG = (^^Q- (35) 

with the auxiliary matrices given by 

u;fzdz + g+{z + dz) -f±fTpa , 
± I _^±^TpG Ufzdz + g4^ + dz' ^ ^ 
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In other words, the eigenvalue problem is reduced to the form 

i^JtiSVxJlSV (37) 





that can be written as the coupled differential equations, 

[Ufzd, + g+{z + d,) - x±] 0^ - y0f ± y 0f = 0, (38) 

[cofzd, + g4^ + d,) - x±] <Pt - ± = 0, (39) 
and their sign inversion counterparts, 

[cofzd, - g+iz + d^) - x±] - Y^t ± = 0, (40) 

[ujzd, -g-{z + d^) - x±] - Y^f ± = 0, (41) 

where we have used the shorthand notation 0^ = (pfiz) and (p^ = (j)^{—z). At this point 
we can take advantage of parity and use a pair of Bargmann functions with well defined 
parity, $^_|_ = (pf ±(j)f such that Tpc^^^^ = ±^^^^, that yield a coupled differential set, 

[cofzd. - X±] + 9-,iz + 9,)$t_ - ^^^^<^t^ = 0, (42) 

[i^fzd. - X±] + g-{^ + 5.)$t - = 0, (43) 

[u^fzd. - X±] $ t- + 9+{z + 9.)$t+ - = 0, (44) 

l^fzd, - X±] +g4^ + 5.)$ J+ - ^^^^^<f = 0, (45) 
that can be solved by Frobenius method with the power series solutions, 

oo 

$t=E4-''' (46) 



k=0 

oo 

^t=E4-.i-''^\ (47) 

k=0 

for one of the pairs of even and odd parity Bargmann functions; the coefficients for the 
other pair, $^_|_, are obtained as functions of the coefficients from (42) and (44) . 



Such an approach yields a five-term recurrence relation for the coefficients ct 



with 



4(2)4 + af{2)cf + 4(2)4 = 0, (48) 
4(3)C3± + 4(3)4 + 4(3)4 + 4(3)4 = 0, (49) 

4(i)cf + 4(j)cf-l + 4(j)cf-2 + 4(i)cf-3 + 4(j)Cj'-4 = 0, j = 4, 5, . . . 

(50) 



«o (^■) = (J3^ 3+9- h T {-iyu^2] , (51) 

afif) = (j - 1) {g^ [ui T {-iyoo2] [{j - Ipf - x±] + 

+g+ [u, ± i-iyu2\ [x± - U - 2)uf] } , (52) 
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afU) = g+ [ur ± {-lyu^] [x± - (j - 2)ujf] - 
- 9- K T (-!)■' t^2] [x± - (j - 3)u;/] , 



(53) 



a 



± 



(j) = 9+9- K T (-1)^W2] 



Notice that for the case of identical qubits, ui 
the if+ block, 

{ufzd, - x+)'^t+ 



UJ2 



(54) 
(55) 

loq ,the differential set for 



0, 
0, 
0, 
0, 



(56) 
(57) 
(58) 
(59) 



immediately shows us that the Bargman function 0^(-2) has well defined odd parity as 



Eq.(59) accepts only the trivial solution $2+(2;) = 0. In this case the five-term recurrence 
relations reduce to three-term recurrence relations for the coefficients, where we have 
rearranged the coefficients for the sake of simplicity, 
1 



4 



J 



[x± - (j - l)w/] 



2 



J 



1,2,... 



^+[x± - (j - l)c^/] 
Negative superindex corresponds to the differential set for the block, 

{ufzd,-X-)'^i,-+9+{z + d,)<i>^^^ =0, 

{ujfzd,-X-)'^2,+ -^o'^i,+ =0, 
{ufzd, - X-)'^l+ + 9+{z + <9.)$7 - wo$2,+ = 0, 



{ufzd, - X-)'^2 



0. 



(60) 
(61) 

(62) 
(63) 
(64) 
(65) 



with = Sjlo^2j-^^"' ^^"^ ^i- ~ Ylf=o(^%+i^'^^^'^ ■ Notice that the coefficients in the 



three-term recurrence relations satisfy lim 



1. 



4. Time Evolution 

We want to include a way of calculating the exact time evolution in the RWA and 
show that numerical solutions for the full quantum problem are simplified by using the 



parity bases introduced in (21) and (22). In order to round up our example we show 



some typical measurements on radiation- matter interaction systems; e.g., mean photon 



number, n, population inversion, von Neumann entropy of the 



.(2)^ 



bipartite system, S = —pqlnpg, and concurrence as defined by Wootters [H 
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Let us start from the Tavis-Cummings Hamiltonian for two-qubits in the frame defined 
by the transformation Ujij(t) = e^**^/^*, 



H 



RWA 



i=i,2 



(66) 



i=i,2 



[UJi 



where the detunnings are given by Aj 
else [B] that the right unitary transformation, 



Uf) /2. We have pointed out somewhere 



T- 




(67) 



diagonahzes the single-qubit Jaynes-Cummings model in the field basis. So, one can use 
that transformation for each qubit and obtain, 

Ti ® T2) Hsc (fl ® ff 



Hrwa, 



(6J 



where 



( Ai + A2 g2^n - 1 gi\^n - 1 

g2^n - 1 Ai - A2 

gi^/n - 1 -A1 + A2 
gi\/n g2y/n 





-A1-A2 / 

In order to obtain the evolution for this system, we can see that 



(69) 



fl ® T2) Hsc (t1 ® f /) (f 1 ® f^) Hsc (t1 ® fl 
fl ® fa) f f j ® fl) - (fl ® fa) HscMsc ffl 



^ft 



Ti ® T2 ) f f J 



with 



f |0)(0| + |1)(1| \ 

|0)(0| 

|0)(0| 

^ 

which leads us to write the evolution operator of the system in the RWA as 



v 



(70) 

,(71) 
(72) 

(73) 



U 



RWA 



Ti ® T2 ) e"*^^^* i f j 



(74) 



This requires just to keep track of the action of ^f j ® fi j over the given initial state 
in order to calculate its evolution. It is simple to calculate the exponential because the 
Hamiltonian Hgc has a depressed quartic as characteristic polynomial. 



Co 

Cl 
C2 



A/l 



Co + ciA + C2A^ + A^, 
Hgl-gD + Aj-Al] [{n ~ iM 
2{glA, + glA2), 
il-2n)igf + gl)-2{Aj + Al), 



9l) 



(75) 

Al - A^] , (76) 
(77) 
(78) 
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with roots, 



Ai,: 



P 



2c2 - 2ci/p 



-^3,4 

and parameters, 
P = 



-p ± V -p^ - 2c2 + 2ci/p 



12co + (c2 - c) 
3c 



-(g + 27cl - 72coC2 + 2c: 



nl/3 



q = V(27cf - 72coC2 + 2ci)2 - 4(12co + cj)^. 



(79) 
(80) 

(81) 
(82) 
(83) 



Once the eigenvalues are obtained, it is trivial but cumbersome to calculate the 
eigenstates of Hsc that allow us to write the evolution operator U(t). For the sake 
of giving an example, let us take the peculiar and simple case when the initial state is 
a Fock state with the two-qubits in the ground state, |V^(0)) = \n,g,g), with symmetric 
detunings, Ai = — A2 = A, and identical couplings, gi = g2 = g, two of the eigenvalues 
are zero and the other two are X± = ±2a/A2 + g'^{n — 1/2). So, we can write the time 



evolution of the state as 



^^|gp(cos^^(. 



2)t-l)V^ 



\ 



gyn+i 

Q2(n+1) 



[A (cosfi(n + l)t - 1) - ^n{n + 1) smn{h + l)t] V 
j^gj^ [-A (cos n{n + l)t - 1) - ^n{n + 1) sin n{n + l)t] V 



\ 



Q2(n) 



[2A2 + g^{h - 1) + g'^n cosQ{h)t] 



J 



with the Rabi frequency defined as Q{n) = \X±\ = 2a/A^ + g'^{h — 1/2). Notice that 
these initial conditions give = {^f{y)) ^"^^ {^^^^) = ~{(^x'^)- 

For more general parameter sets in the weak-coupling regime, let us consider an 
initial state example given by a coherent field plus the qubits in the ground state, 
I'lpiO)) = \a,g,g), the transformation applied over this state acts like the identity, 
f^^fl) IV^(O)) = IV^(O)). In Fig. 



we show the time evolution for the mean photon 
number, population inversion, von^Nfeuman entropy and qubit bipartite concurrence 
for the case of two nonidentical qubits, ui = l.lw/ and 0J2 = O.Suf, weakly coupled 
to the field, gi = O.OSuf and g2 = O.OAuf. One can see that the collapse and revival 
of the population inversion, characteristic of such an initial state, is not as strong as 
in the case of equal qubits [H]. The fast modulation of the curves is an interesting 
feature that does not occur in the identical qubits case [H]. The von Neuman entropy 
points to entanglement between the qubit ensemble and the quantum field. The qubits 
bipartite concurrence tells us that quantum correlations between the qubits are created 
by interaction with the field but the qubit ensemble is close to being separable. 
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t (units of hiOf 



1000 



Figure 4. (Color online) The time evolution of the (a) mean photon number, (b) 
population inversion, (c) von Neuman entropy and (d) bipartite concurrence for a 
quantum Rabi Hamiltonian under the RWA with the parameter set {1^1,^2,51,52} = 
{1.1, 0.3, 0.03, 0.04}w/ and initial state |0(O)) = \a,g,g) with a = V^. 



4-2. Parity based numeric approach 

In the most general case of any given coupling parameter set, one can approximate 



results by adequately truncating the H± matrices in (23) obtained by using the parity 



bases (21 ) and (22 ). This numerical approach allows us to calculate the time evolution of 



a general initial state decomposed in the even and odd parity chain bases, respectively. 



n=0 

00 

n=0 



.(+) 



.(+) 
-'in 



W0))-= 



n=0 

00 

n=0 



.(-) 



.(-) 
^4n 



"^n^g^g) + ci;^^i|2n, e, e) + 

2n+ l,e,^) +4t|3|2n+ l,^,e), 

n)-, 

2n,e,g) + ci~l-^\2n, g,e) + 

2n + 1, ^, g) + c^3|2n + 1, e, e) 



with the prescription, 

c-<^)(t) = 



c-^"^(0). 



^5) 



^6) 
^7) 



^9) 
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where the notation c'-^^ corresponds to a vector of state amphtudes. The average photon 
number of the field is then given by 



im) = {Ht))+ + (nit))- 



.(±)|2 
^4n I 



+ C, 



(i) I 
4n+ll 



n=0 



+ {2n + l) (^1 



,(±) I 

4n+2 I 



+ c 



(±) I 
4n+3 I 



(90) 

.(91) 



any other quantity of interest can be calculated from the qubit ensemble reduced density 



matrix, = 



f{n\p\n)f 
( 



n=0 



E 

n=0 



C 



(+) 



^4ra+l I 
(+) (+)* 
'^4n+2'^4n+l 
(+) (+)* 

J+)J+)* 
"^4n ^An+l 

!„(-) 12 
I '-4n+3 1 

„(-)„(-)* 

(-)* (-) 

(-) J-)* 
4n+2'-4n+3 



with p 

„(+)* „(+) 

|J+) 12 
I '-4n+2 I 

^4n+3^4:n+2 

'-4n ^4n+2 

„(-) „(-)* 
^4n+3^4n 

|J-)|2 
l"^4n I 

"^4n+l"^4n 

J-)* 
'^4n+2'^4n 



in the a^z'^ (g) 0"^^'' basis, 



(+) ^(+) ^(+)* 



''4ra+l'^4n+3 

^(+) ^(+)* 
'-4n+2'-4n+3 

12 

I '-4?i+3 1 

<-4n '-4n+3 

(-) {-)* 
^4n+3^4n+l 

„(-)„(-)* 
"^4n "^4n+l 



C 



(-) 
4n+l I 



(-) {-)* 
'^4n+2'^4ri+l 



'^4n+1^4n 

S+) 1+)* 

^4n+3'-'4n 
|J+)|2 
\^4n I 

'^4n+3^4n+2 

„(-)„(-)* 
'-4n ^4n+2 

.(-) J-)* 
.(-) 



(92) 



'^4n+l'^4?i+2 



(93) 



c 



4n+2 1 



In Fig. |5] we show a comparison between the dynamics given by the RWA and the full 
quantum Rabi Hamiltonian for coupling parameters near to 10% of the field frequency 
value. One can see that the RWA fails in providing a detailed description of the dynamics 
as mentioned above and in Figure |6] presents the dynamics given by the quantum 
Rabi Hamiltonian for strong-coupling parameters, one can see a great deviation in the 
behavior from the weak-coupling dynamics presented in Fig. |4] where the couplings are 
a hundred times smaller. Note that the quantum correlations between the qubit and 
the field and between the qubits increase with respect to the weak-couphng case. 



5. Conclusion 



This manuscript was motivated by the results of [3l], where the assumption 1^2 = is 
made in order to deliver three-term recurrence relations and Green functions similar in 
form to those presented by [26] for the single-qubit case. We wanted to show that it is 
possible to study a proper, physical two-qubit quantum Rabi Hamiltonian by extending 
well known methods for single-qubit models that exist in the literature since [HIIS]- 

We have shown the well-known result that the spectra of the two-qubit quantum 
Rabi model in the weak-coupling regime is easily obtained by partitioning the Hilbert 
space in manifolds of constant total number of excitation, a standard procedure in 
the single-qubit quantum Rabi problem [10]. Our main results are the following. 
In the strong-coupling regime we have calculated the eigenvalues up to second order 
perturbation correction. An interesting result occurs for ultrastrong-couplings where 
the eigenvalues of the two-qubit Rabi model are described by two forced oscillator 
spectral series; this is equivalent to what happens in the single-qubit case [HIIH]- In the 
intermediate regime, partitioning the Hilbert space in even and odd parity subspaces 
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t (units of hojf) t (units of huj) 



Figure 5. (Color online) The time evolution of the (a) mean photon number, (b) 
population inversion, (c) von Neuman entropy and (d) bipartite concurrence for a 
RWA Hamiltonian (dashed blue) and a quantum Rabi Hamiltonian (solid red) with the 
parameter set {uji,uj2, 91,92} = {1-1, 0.3, 0.1, O.lSjwj and initial state 10(0)) = \a,9,9) 
with a = \/2. 

simplifies the numerical calculation of the eigenvalues; an interesting result from numeric 
evidence is that the ground state of even and odd subspaces is degenerate for coupling 
values larger than the field frequency value. The numerical spectra points to energy 
crossings within and between each of the parity subspaces with moderate couplings; 
the case studied in [Ml seems to be an exception. The crossings within the spectra 
of parity subspaces represent a different behavior to that of the single-qubit results, 
where the spectra branches within each of the parity subspaces do not cross and the 
parity subspaces spectra cross each other [111 126]. We have obtained the normal 
modes for the two-qubit Rabi model by linear algebra methods via parity subspaces 
as four-terms recurrence relations for the coefficients of the eigenstates; in the case of 
identical couplings the linear system becomes singular and we cannot provide an analytic 
solution, nevertheless it is possible to obtain numeric solutions for this particular case. In 
Bargmann representation, we find that the recurrence relation for the continued proper 
functions coefficients are five-term and reduce to three-term recurrence relations in the 
case of identical qubits; additionally, for identical qubits, one of the proper functions has 
well defined parity. Finally, we wanted to show how powerful the parity decomposition 
is by calculating the time evolution of quantities related to the degree of entanglement 
in the ensemble-field and qubit system for a simple initial state of the qubits in their 
ground state and the field in a coherent state with two photons on average under weak-. 
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t (units of hojf) t (units of huj) 

Figure 6. (Color online) The time evolution of the (a) mean photon number, (b) 
population inversion, (c) von Neuman entropy and (d) bipartite concurrence for a 
quantum Rabi Hamiltonian with the parameter set {oji,uj2, 9i, 92} = {1-1, 0.3, 3, 4}wj 
and initial state \4'{0)) — \a,g,g) with a — \/2 considering a subspace of up to one 
thousand excitations. 

moderate- and strong-coupling. 
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